# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools", "sdmTMB", "sdmTMBextra", "terra", "mapplots",
"viridis", "visreg", "modelr", "future", "kableExtra", "ggh4x", "patchwork")
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
}
invisible(lapply(pkgs, library, character.only = T))
# Import some plotting functions
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")
# devtools::install_github("seananderson/ggsidekick") # not on CRAN
library(ggsidekick)
theme_set(theme_sleek())
# Set path
home <- here::here()
# For crossvalidation: paralell processing
plan(multisession)
set.seed(99) Fit stomach content models
Load libraries
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/main-analysis/03-fit-diet-models_cache/html"))Read cleaned data
df <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |>
#filter(pred_length > 15 & pred_length < 50) |>
mutate(depth_sc = (depth - mean(depth))/sd(depth),
year_f = as.factor(year),
month_f = as.factor(month),
ices_rect = as.factor(ices_rect),
pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length)) # read_csv(paste0(home, "/data/clean/stomachs.csv")) |> summarise(mean = mean(pred_length))
# Read the prediction grid...
pred_grid <- bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))
# Summarize samples size per location for plotting...
dd_plot <- df |>
group_by(year, X, Y) |>
mutate(sample_size = n(),
pos_id = paste(year, X, Y)) |>
ungroup() |>
distinct(pos_id, .keep_all = TRUE)
# Scale with respect to data!
df |> group_by(month) |> summarise(n = n()) |> arrange(desc(n))# A tibble: 8 × 2
month n
<dbl> <int>
1 3 7040
2 11 2207
3 12 1980
4 4 553
5 2 549
6 10 57
7 6 36
8 5 3
pred_grid <- pred_grid |>
drop_na(depth) |>
filter(depth < quantile(df$depth, probs = 0.99)) |>
filter(quarter == 4) |> # Not needed in theory for saduria...
mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
year_f = as.factor(year),
month_f = as.factor(11),
year = as.integer(year),
pred_length_sc = 0,
ices_rect = as.factor(ices_rect)) |>
droplevels()
ggplot() +
geom_histogram(data = pred_grid, aes(depth, fill = "pred_grid")) +
geom_histogram(data = df, aes(depth, fill = "data"))ggplot(pred_grid |> filter(year == 1999), aes(X, Y, color = depth)) +
geom_point()Plot data!
df |>
pivot_longer(c(FR_spr, FR_her, FR_sad)) |>
ggplot(aes(year, value)) +
geom_jitter(height = 0, alpha = 0.5) +
coord_cartesian(ylim = c(0, 0.1)) +
facet_wrap(~name, ncol = 1) +
geom_smooth()Set up a mesh
mesh <- make_mesh(df, c("X", "Y"), cutoff = 6)
ggplot() +
inlabru::gg(mesh$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = df, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) +
labs(x = "Easting (km)", y = "Northing (km)")Spatial cross validation and AIC to compare spatial vs non spatial models
5-fold spatial cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects
# Set up spatial clusters
clust <- kmeans(df[, c("X", "Y")], 5)$cluster
# Plot
df$clust_id <- clust
plot_map +
geom_point(data = df, aes(X*1000, Y*1000, color = as.factor(clust_id))) +
scale_color_viridis(discrete = TRUE) +
geom_sf(size = 0.1)Warning: Removed 39 rows containing missing values (`geom_point()`).
Warning: colourbar guide needs continuous scales.
# Sprat
# ICES as a random effects
spr_space_1_cv <- sdmTMB_cv(
FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
spr_space_2_cv <- sdmTMB_cv(
FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Herring
# ICES as a random effects
her_space_1_cv <- sdmTMB_cv(
FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
her_space_2_cv <- sdmTMB_cv(
FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Saduria
# ICES as a random effects
sad_space_1_cv <- sdmTMB_cv(
FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
spatial = "off",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
# Spatial random field
sad_space_2_cv <- sdmTMB_cv(
FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
spatial = "on",
family = tweedie(link = "log"),
fold_ids = clust,
k_folds = length(unique(clust))
)Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.
AIC
# ICES as a random effects
fit_spr_m1 <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_her_m1 <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
fit_sad_m1 <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "off",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_sad_m1)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_smooth_sigma` standard error may be large
ℹ Try simplifying the model, adjusting the mesh, or adding priors
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
# Spatial random field
fit_spr_m2 <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_spr_m2)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
fit_her_m2 <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
sanity(fit_her_m2)✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
fit_sad_m2 <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
Print AIC & cross validation table
# Higher is better!
# likelihood from spatial cv
sum_loglik_space <- data.frame(Prey = c("Sprat", "Herring", "Saduria"),
rec = c(spr_space_1_cv$sum_loglik / 10000000,
her_space_1_cv$sum_loglik / 10000000,
sad_space_1_cv$sum_loglik / 10000000),
space = c(spr_space_2_cv$sum_loglik / 10000000,
her_space_2_cv$sum_loglik / 10000000,
sad_space_2_cv$sum_loglik / 10000000)) |>
mutate(rec_temp = rec, space_temp = space,
rec = ifelse(rec_temp > space_temp, paste0("**", rec, "**"), rec),
space = ifelse(space_temp > rec_temp, paste0("**", space, "**"), space)) |>
rename("sum loglik<sub>ICES rect</sub>" = rec,
"sum loglik<sub>spatial</sub>" = space) |>
dplyr::select(-rec_temp, -space_temp)
sum_loglik_space Prey sum loglik<sub>ICES rect</sub> sum loglik<sub>spatial</sub>
1 Sprat **5.58014982903663e-05** 4.90488173952715e-05
2 Herring -0.850543167016049 **-0.807858283667988**
3 Saduria -13.6423530384996 **-9.92716500451823**
# AIC ices models
aic_ices <- AIC(fit_spr_m1, fit_her_m1, fit_sad_m1) |>
rownames_to_column() |>
mutate(Prey = "Sprat",
Prey = ifelse(rowname == "fit_her_m1", "Herring", Prey),
Prey = ifelse(rowname == "fit_sad_m1", "Saduria", Prey)) |>
mutate(AIC_r = AIC,
AIC_r2 = AIC) |>
dplyr::select(-rowname, -AIC, -df)
aic_space <- AIC(fit_spr_m2, fit_her_m2, fit_sad_m2) |>
rownames_to_column() |>
mutate(Prey = "Sprat",
Prey = ifelse(rowname == "fit_her_m2", "Herring", Prey),
Prey = ifelse(rowname == "fit_sad_m2", "Saduria", Prey)) |>
mutate(AIC_s = AIC,
AIC_s2 = AIC) |>
dplyr::select(-rowname, -AIC, -df)
aic <- aic_ices |>
left_join(aic_space, by = "Prey") |>
mutate(AIC_r = ifelse(AIC_r2 < AIC_s2, paste0("**", AIC_r, "**"), AIC_r),
AIC_s = ifelse(AIC_s2 < AIC_r2, paste0("**", AIC_s, "**"), AIC_s)) |>
rename("AIC<sub>spatial</sub>" = AIC_s,
"AIC<sub>ICES rect</sub>" = AIC_r) |>
dplyr::select(-AIC_s2, -AIC_r2)
aic Prey AIC<sub>ICES rect</sub> AIC<sub>spatial</sub>
1 Sprat -2925.7767 **-3284.89068021079**
2 Herring 558.8468 **426.356157207818**
3 Saduria -1323.1981 **-1547.59615868862**
kableExtra::kbl(sum_loglik_space, format = "markdown", caption = "5-fold Spatial cross validation", digits = 2) |>
kable_styling(font_size = 20)| Prey | sum loglikICES rect | sum loglikspatial |
|---|---|---|
| Sprat | 5.58014982903663e-05 | 4.90488173952715e-05 |
| Herring | -0.850543167016049 | -0.807858283667988 |
| Saduria | -13.6423530384996 | -9.92716500451823 |
kableExtra::kbl(aic, format = "markdown", caption = "AIC", digits = 2) |>
kable_styling(font_size = 20)| Prey | AICICES rect | AICspatial |
|---|---|---|
| Sprat | -2925.78 | -3284.89068021079 |
| Herring | 558.85 | 426.356157207818 |
| Saduria | -1323.20 | -1547.59615868862 |
# TODO: it seems perhaps ?kbl doesn't obey rounding when I use **. Make table in word instead for the reportCheck residuals and from selected models
# Summary
summary(fit_spr_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_spr ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc 0.37 0.07
spred_length_sc 0.52 0.07
Smooth terms:
Std. Dev.
sds(pred_length_sc) 35.52
Random intercepts:
Std. Dev.
month_f 0.82
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -6.45 0.43
(Intercept)-1994 -6.22 0.46
(Intercept)-1995 -7.38 0.43
(Intercept)-1996 -6.08 0.39
(Intercept)-1997 -7.31 0.45
(Intercept)-1998 -5.85 0.49
(Intercept)-1999 -6.20 0.38
(Intercept)-2000 -6.27 0.39
(Intercept)-2001 -5.66 0.40
(Intercept)-2002 -7.70 0.46
(Intercept)-2003 -4.69 0.40
(Intercept)-2004 -7.35 0.39
(Intercept)-2005 -6.20 0.38
(Intercept)-2006 -6.14 0.38
(Intercept)-2007 -6.35 0.37
(Intercept)-2008 -6.25 0.37
(Intercept)-2009 -5.90 0.38
(Intercept)-2012 -4.94 0.38
(Intercept)-2013 -5.62 0.37
(Intercept)-2014 -5.16 0.42
(Intercept)-2018 -6.00 0.39
(Intercept)-2021 -5.88 0.37
Dispersion parameter: 0.30
Tweedie p: 1.44
Matérn range: 16.14
Spatial SD: 1.24
ML criterion at convergence: -1651.445
See ?tidy.sdmTMB to extract these values as a data frame.
summary(fit_her_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_her ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc -0.23 0.09
spred_length_sc -0.79 0.07
Smooth terms:
Std. Dev.
sds(pred_length_sc) 34.7
Random intercepts:
Std. Dev.
month_f 0.12
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -9.63 0.49
(Intercept)-1994 -8.07 0.45
(Intercept)-1995 -7.29 0.37
(Intercept)-1996 -8.16 0.37
(Intercept)-1997 -6.92 0.35
(Intercept)-1998 -6.93 0.46
(Intercept)-1999 -7.62 0.31
(Intercept)-2000 -7.49 0.33
(Intercept)-2001 -7.40 0.35
(Intercept)-2002 -8.26 0.47
(Intercept)-2003 -5.66 0.28
(Intercept)-2004 -7.27 0.28
(Intercept)-2005 -8.15 0.31
(Intercept)-2006 -7.49 0.27
(Intercept)-2007 -7.77 0.28
(Intercept)-2008 -7.48 0.26
(Intercept)-2009 -6.65 0.29
(Intercept)-2012 -6.35 0.32
(Intercept)-2013 -6.29 0.28
(Intercept)-2014 -7.68 0.49
(Intercept)-2018 -8.13 0.32
(Intercept)-2021 -8.17 0.25
Dispersion parameter: 0.61
Tweedie p: 1.48
Matérn range: 16.11
Spatial SD: 1.28
ML criterion at convergence: 204.178
See ?tidy.sdmTMB to extract these values as a data frame.
summary(fit_sad_m2)Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc -0.01 0.15
spred_length_sc 0.19 0.05
Smooth terms:
Std. Dev.
sds(pred_length_sc) 0
Random intercepts:
Std. Dev.
month_f 0.95
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -9.19 1.08
(Intercept)-1994 -9.22 1.08
(Intercept)-1995 -8.63 1.01
(Intercept)-1996 -8.61 1.00
(Intercept)-1997 -8.61 1.01
(Intercept)-1998 -9.71 1.11
(Intercept)-1999 -8.58 0.99
(Intercept)-2000 -8.16 0.98
(Intercept)-2001 -7.53 0.98
(Intercept)-2002 -9.07 1.00
(Intercept)-2003 -8.43 0.98
(Intercept)-2004 -8.44 0.98
(Intercept)-2005 -8.24 0.98
(Intercept)-2006 -8.01 0.97
(Intercept)-2007 -8.80 0.98
(Intercept)-2008 -8.67 0.98
(Intercept)-2009 -8.61 0.99
(Intercept)-2012 -7.96 0.99
(Intercept)-2013 -8.48 0.99
(Intercept)-2014 -8.22 0.98
(Intercept)-2018 -7.54 0.95
(Intercept)-2021 -7.14 0.94
Dispersion parameter: 0.73
Tweedie p: 1.54
Matérn range: 80.79
Spatial SD: 2.99
ML criterion at convergence: -782.798
See ?tidy.sdmTMB to extract these values as a data frame.
**Possible issues detected! Check output of sanity().**
# Residuals
spr_res <- mcmc_res <- residuals(fit_spr_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_spr_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.002765 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.65 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 32.4323 seconds (Warm-up)
Chain 1: 0.064402 seconds (Sampling)
Chain 1: 32.4967 seconds (Total)
Chain 1:
df$spr_res <- as.vector(spr_res)
her_res <- mcmc_res <- residuals(fit_her_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_her_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.002747 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.47 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 19.0425 seconds (Warm-up)
Chain 1: 0.030647 seconds (Sampling)
Chain 1: 19.0732 seconds (Total)
Chain 1:
df$her_res <- as.vector(her_res)
sad_res <- mcmc_res <- residuals(fit_sad_m2, type = "mle-mcmc",
mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_sad_m2,
mcmc_iter = 201,
mcmc_warmup = 200))
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00279 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.9 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 201 [ 0%] (Warmup)
Chain 1: Iteration: 20 / 201 [ 9%] (Warmup)
Chain 1: Iteration: 40 / 201 [ 19%] (Warmup)
Chain 1: Iteration: 60 / 201 [ 29%] (Warmup)
Chain 1: Iteration: 80 / 201 [ 39%] (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%] (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%] (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%] (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%] (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%] (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%] (Warmup)
Chain 1: Iteration: 201 / 201 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 243.25 seconds (Warm-up)
Chain 1: 0.251266 seconds (Sampling)
Chain 1: 243.501 seconds (Total)
Chain 1:
df$sad_res <- as.vector(sad_res)
# Plot all together
df |>
rename(Sprat = spr_res,
Herring = her_res,
Saduria = sad_res) |>
pivot_longer(c(Sprat, Herring, Saduria)) |>
ggplot(aes(sample = value)) +
facet_wrap(~name) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)rename: renamed 3 variables (Sprat, Herring, Saduria)
pivot_longer: reorganized (Sprat, Herring, Saduria) into (name, value) [was 12425x25, now 37275x24]
ggsave(paste0(home, "/figures/supp/diet_qq.pdf"), width = 17, height = 11, units = "cm")Plot conditional effects, random effects and spatial predictions
# Depth
vis_spr_dep <- visreg(fit_spr_m2, xvar = "depth_sc", plot = FALSE)
vis_her_dep <- visreg(fit_her_m2, xvar = "depth_sc", plot = FALSE)
vis_sad_dep <- visreg(fit_sad_m2, xvar = "depth_sc", plot = FALSE)
vis_dep <- bind_rows(vis_spr_dep$fit |> mutate(species = "Sprat"),
vis_her_dep$fit |> mutate(species = "Herring"),
vis_sad_dep$fit |> mutate(species = "Saduria")) |>
mutate(var = "Depth (scaled)") |>
rename(x = depth_sc)
d_dep <- bind_rows(vis_spr_dep$res |> mutate(species = "Sprat"),
vis_her_dep$res |> mutate(species = "Herring"),
vis_sad_dep$res |> mutate(species = "Saduria")) |>
mutate(var = "Depth (scaled)") |>
rename(x = depth_sc)
# Month
vis_spr_mon <- visreg(fit_spr_m2, xvar = "month_f", plot = FALSE)
vis_her_mon <- visreg(fit_her_m2, xvar = "month_f", plot = FALSE)
vis_sad_mon <- visreg(fit_sad_m2, xvar = "month_f", plot = FALSE)
vis_mon <- bind_rows(vis_spr_mon$fit |> mutate(species = "Sprat"),
vis_her_mon$fit |> mutate(species = "Herring"),
vis_sad_mon$fit |> mutate(species = "Saduria")) |>
mutate(var = "Month") |>
rename(x = month_f) |>
mutate(x = as.numeric(as.character(x)))
d_mon <- bind_rows(vis_spr_mon$res |> mutate(species = "Sprat"),
vis_her_mon$res |> mutate(species = "Herring"),
vis_sad_mon$res |> mutate(species = "Saduria")) |>
mutate(var = "Month") |>
rename(x = month_f) |>
mutate(x = as.numeric(as.character(x)))
# Predator length
vis_spr_len <- visreg(fit_spr_m2, xvar = "pred_length_sc", plot = FALSE)
vis_her_len <- visreg(fit_her_m2, xvar = "pred_length_sc", plot = FALSE)
vis_sad_len <- visreg(fit_sad_m2, xvar = "pred_length_sc", plot = FALSE)
vis_len <- bind_rows(vis_spr_len$fit |> mutate(species = "Sprat"),
vis_her_len$fit |> mutate(species = "Herring"),
vis_sad_len$fit |> mutate(species = "Saduria")) |>
mutate(var = "Predator length") |>
rename(x = pred_length_sc)
d_len <- bind_rows(vis_spr_len$res |> mutate(species = "Sprat"),
vis_her_len$res |> mutate(species = "Herring"),
vis_sad_len$res |> mutate(species = "Saduria")) |>
mutate(var = "Predator length") |>
rename(x = pred_length_sc)
vis <- bind_rows(vis_dep, vis_mon, vis_len)
vis_dat <- bind_rows(d_dep, d_mon, d_len)
# Plot!
ggplot(vis |> filter(!var == "Month"), aes(x = x, y = visregFit)) +
facet_grid(var ~ species) +
geom_point(data = vis_dat |> filter(!var == "Month"), aes(x = x, y = visregRes),
alpha = 0.3, color = "gray50", size = 0.2) +
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.3, color = NA) +
geom_line(color = "steelblue", linewidth = 1) +
labs(x = "Scaled variable", y = "Prediction") +
NULLggsave(paste0(home, "/figures/supp/conditional.pdf"), width = 17, height = 11, units = "cm")
# Now do month (categorial)
ggplot(vis |> filter(var == "Month"), aes(x = as.factor(x), y = visregFit)) +
facet_wrap(~ species) + # free grid?
geom_jitter(data = vis_dat |> filter(var == "Month"), aes(x = as.factor(x), y = visregRes),
alpha = 0.15, color = "gray60", size = 0.2, height = 0) +
geom_errorbar(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.6, color = "steelblue",
width = 0, linewidth = 1) +
geom_point(color = "steelblue", size = 2) +
labs(x = "Month", y = "Prediction") +
theme(aspect.ratio = 5/6) +
NULLggsave(paste0(home, "/figures/supp/conditional_month.pdf"), width = 17, height = 6, units = "cm")Calculate indices
# Need to refit the models with the extra time argument. The reason we don't do that before is because it creates a mismatch in residual dimensions and data
# This is for interpolating between year using the random walk
extra_time <- pred_grid |> filter(!year %in% df$year) |> distinct(year) |> pull(year)filter: removed 322,828 rows (76%), 102,718 rows remaining
distinct: removed 102,711 rows (>99%), 7 rows remaining
# Spatial random field
fit_spr_m2_extra <- sdmTMB(
FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
fit_her_m2_extra <- sdmTMB(
FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
fit_sad_m2_extra <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
extra_time = extra_time,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))
# Predict on grid, for indices and maps
pred_spr <- predict(fit_spr_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_her <- predict(fit_her_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_sad <- predict(fit_sad_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
# Make temporal index!
ncells <- filter(pred_grid, year == max(pred_grid$year)) |> nrow()filter: removed 410,872 rows (97%), 14,674 rows remaining
index_spr <- get_index(pred_spr, area = rep(1/ncells, nrow(pred_spr$data)), bias_correct = TRUE)
index_her <- get_index(pred_her, area = rep(1/ncells, nrow(pred_her$data)), bias_correct = TRUE)
index_sad <- get_index(pred_sad, area = rep(1/ncells, nrow(pred_sad$data)), bias_correct = TRUE)
# Make long
index <- bind_rows(index_spr |> mutate(prey = "Sprat"),
index_her |> mutate(prey = "Herring"),
index_sad |> mutate(prey = "Saduria")) |>
mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA
# For comparison with data
df_sum <- df |>
pivot_longer(c("FR_spr", "FR_her", "FR_sad"), names_to = "prey", values_to = "fr") |> # because the model cannot predict the large tails!
group_by(prey) |>
mutate(upr = quantile(fr, probs = 0.99)) |>
ungroup() |>
filter(fr < upr) |>
group_by(year, prey) |>
summarise(mean_fr = mean(fr, na.rm = TRUE)) |>
mutate(prey = ifelse(prey == "FR_spr", "Sprat", prey),
prey = ifelse(prey == "FR_her", "Herring", prey),
prey = ifelse(prey == "FR_sad", "Saduria", prey))pivot_longer: reorganized (FR_sad, FR_spr, FR_her) into (prey, fr) [was 12425x25, now 37275x24]
group_by: one grouping variable (prey)
mutate (grouped): new variable 'upr' (double) with 3 unique values and 0% NA
ungroup: no grouping variables
filter: removed 375 rows (1%), 36,900 rows remaining
group_by: 2 grouping variables (year, prey)
summarise: now 66 rows and 3 columns, one group variable remaining (year)
mutate (grouped): changed 66 values (100%) of 'prey' (0 new NA)
index |>
ggplot(aes(year, est, fill = Observed)) +
geom_point(shape = 21, alpha = 0.7) +
scale_fill_manual(values = c("white", "grey10")) +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
labs(x = "Year", y = "Per capita predation", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme_sleek(base_size = 9) +
theme(legend.position = "bottom")ggsave(paste0(home, "/figures/supp/index_ci.pdf"), width = 17, height = 7, units = "cm")
index |>
mutate(upr = ifelse(prey == "Sprat" & upr > 0.015, 0.015, upr)) |>
ggplot(aes(year, est)) +
geom_point(alpha = 0.7) +
stat_smooth(method = "gam", formula = y~s(x, k=3), color = "steelblue") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
labs(x = "Year", y = "Per capita predation", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))mutate: changed 7 values (8%) of 'upr' (0 new NA)
ggsave(paste0(home, "/figures/index.pdf"), width = 17, height = 7, units = "cm")Sensitivity with respect to prediction grid
nd_year <- data.frame(year = unique(pred_grid$year),
depth_sc = 0,
pred_length_sc = 0,
month_f = as.factor(11))
year_pred_spr <- predict(fit_spr_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)
year_pred_her <- predict(fit_her_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)
year_pred_sad <- predict(fit_sad_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)
year_cond <- bind_rows(year_pred_spr |> mutate(prey = "Sprat"),
year_pred_her |> mutate(prey = "Herring"),
year_pred_sad |> mutate(prey = "Saduria"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
year_cond <- year_cond |> mutate(data = ifelse(year %in% extra_time, "Missing", "Present"))mutate: new variable 'data' (character) with 2 unique values and 0% NA
# Save conditional effect of year
p1 <- ggplot(year_cond |> filter(prey == "Sprat"), aes(year, exp(est), shape = data)) +
geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
color = "steelblue3") +
geom_point(color = "steelblue3") +
#facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) +
labs(x = "", y = "", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
scale_shape_manual(values = c(21, 16)) +
guides(shape = "none") +
coord_cartesian(ylim = c(0, 0.017)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))filter: removed 58 rows (67%), 29 rows remaining
p1p2 <- ggplot(year_cond |> filter(prey == "Herring"), aes(year, exp(est), shape = data)) +
geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
color = "steelblue3") +
geom_point(color = "steelblue3") +
#facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) +
labs(x = "", y = "Per capita predation", shape = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
scale_shape_manual(values = c(21, 16)) +
coord_cartesian(ylim = c(0, 0.008)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))filter: removed 58 rows (67%), 29 rows remaining
p2p3 <- ggplot(year_cond |> filter(prey == "Saduria"), aes(year, exp(est), shape = data)) +
geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
color = "steelblue3") +
geom_point(color = "steelblue3") +
#facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) +
labs(x = "Year", y = "", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
scale_shape_manual(values = c(21, 16)) +
guides(shape = "none") +
coord_cartesian(ylim = c(0, 0.0012)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))filter: removed 58 rows (67%), 29 rows remaining
p3# p1 + p2 + p3
#
# ggsave(paste0(home, "/figures/conditional_year.pdf"), width = 17, height = 7, units = "cm")
p1 / p2 / p3ggsave(paste0(home, "/figures/conditional_year2.pdf"), width = 11, height = 17, units = "cm")
# Predict on full grid (not depth-trimmed) and calculate index
pred_grid2 <-
bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv"))) |>
filter(quarter == 4) |> # Not needed in theory for saduria...
mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
year_f = as.factor(year),
month_f = as.factor(3),
year = as.integer(year),
pred_length_sc = 0,
ices_rect = as.factor(ices_rect)) |>
droplevels()Rows: 457436 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): substrate, month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 490110 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): substrate, month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 473,773 rows (50%), 473,773 rows remaining
mutate: converted 'year' from double to integer (0 new NA)
converted 'ices_rect' from character to factor (0 new NA)
new variable 'depth_sc' (double) with 7,602 unique values and 0% NA
new variable 'year_f' (factor) with 29 unique values and 0% NA
new variable 'month_f' (factor) with one unique value and 0% NA
new variable 'pred_length_sc' (double) with one unique value and 0% NA
# Predict on grid, for indices and maps
pred_spr2 <- predict(fit_spr_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)
pred_her2 <- predict(fit_her_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)
pred_sad2 <- predict(fit_sad_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)
# Make temporal index!
ncells2 <- filter(pred_grid2, year == max(pred_grid2$year)) |> nrow()filter: removed 457,436 rows (97%), 16,337 rows remaining
index_spr2 <- get_index(pred_spr2, area = rep(1/ncells2, nrow(pred_spr2$data)), bias_correct = TRUE)
index_her2 <- get_index(pred_her2, area = rep(1/ncells2, nrow(pred_her2$data)), bias_correct = TRUE)
index_sad2 <- get_index(pred_sad2, area = rep(1/ncells2, nrow(pred_sad2$data)), bias_correct = TRUE)
# Make long
index2 <- bind_rows(index_spr2 |> mutate(prey = "Sprat"),
index_her2 |> mutate(prey = "Herring"),
index_sad2 |> mutate(prey = "Saduria")) |>
mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA
# Calculate the index for the most common ices rectangles, the ones that cumulatively make up most stomach samples!
common_rects <- df |>
mutate(n = n()) |>
group_by(ices_rect) |>
summarise(n_rec = n()) |>
ungroup() |>
arrange(desc(n_rec)) |>
mutate(cumsum = cumsum(n_rec),
sum = sum(n_rec),
perc = cumsum / sum) |>
filter(perc < 0.9)mutate: new variable 'n' (integer) with one unique value and 0% NA
group_by: one grouping variable (ices_rect)
summarise: now 40 rows and 2 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'cumsum' (integer) with 40 unique values and 0% NA
new variable 'sum' (integer) with one unique value and 0% NA
new variable 'perc' (double) with 40 unique values and 0% NA
filter: removed 27 rows (68%), 13 rows remaining
common_rects# A tibble: 13 × 5
ices_rect n_rec cumsum sum perc
<fct> <int> <int> <int> <dbl>
1 42H0 3017 3017 12425 0.243
2 40H0 1573 4590 12425 0.369
3 43H0 1398 5988 12425 0.482
4 43H1 1272 7260 12425 0.584
5 41G9 1211 8471 12425 0.682
6 41H0 630 9101 12425 0.732
7 44H1 605 9706 12425 0.781
8 38G5 317 10023 12425 0.807
9 40G9 315 10338 12425 0.832
10 40G6 245 10583 12425 0.852
11 41G8 220 10803 12425 0.869
12 39G7 178 10981 12425 0.884
13 39G5 146 11127 12425 0.896
pred_grid3 <- filter(pred_grid, ices_rect %in% common_rects$ices_rect)filter: removed 308,038 rows (72%), 117,508 rows remaining
plot_map +
geom_raster(data = pred_grid3 |> filter(year == 2000),
aes(X*1000, Y*1000, fill = ices_rect)) +
geom_point(data = dd_plot, aes(X*1000, Y*1000, size = sample_size), alpha = 0.5) +
scale_size(range = c(.01, 2), name = "# stomachs")filter: removed 113,456 rows (97%), 4,052 rows remaining
Warning: Removed 2 rows containing missing values (`geom_point()`).
Warning: colourbar guide needs continuous scales.
pred_spr3 <- predict(fit_spr_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)
pred_her3 <- predict(fit_her_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)
pred_sad3 <- predict(fit_sad_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)
# Make temporal index!
ncells3 <- filter(pred_grid3, year == max(pred_grid3$year)) |> nrow()filter: removed 113,456 rows (97%), 4,052 rows remaining
index_spr3 <- get_index(pred_spr3, area = rep(1/ncells3, nrow(pred_spr3$data)), bias_correct = TRUE)
index_her3 <- get_index(pred_her3, area = rep(1/ncells3, nrow(pred_her3$data)), bias_correct = TRUE)
index_sad3 <- get_index(pred_sad3, area = rep(1/ncells3, nrow(pred_sad3$data)), bias_correct = TRUE)
# Make long
index3 <- bind_rows(index_spr3 |> mutate(prey = "Sprat"),
index_her3 |> mutate(prey = "Herring"),
index_sad3 |> mutate(prey = "Saduria")) |>
mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA
index |>
ggplot(aes(year, est, fill = Observed)) +
geom_point(shape = 21, alpha = 0.7) +
scale_fill_manual(values = c("white", "grey10")) +
facet_wrap(~prey, ncol = 1, scales = "free") +
geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = year_cond, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_ribbon(data = year_cond, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
geom_point(data = index2, aes(year, est, color = "Full grid index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = index3, aes(year, est, color = "Common grid index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme_sleek(base_size = 9) +
theme(legend.position = "bottom")ggsave(paste0(home, "/figures/supp/index_ci_sens.pdf"), width = 15, height = 30, units = "cm")
# Is the mismatch between data and conditional simply because data happen to be in saduria-rich areas?
saduria <- terra::rast(paste0(home, "/data/saduria-tif/FWBiomassm_raster_19812019presHighweightcor_no0_newZi.tif"))
WGS84 <- "+proj=longlat +datum=WGS84"
saduria_latlon <- terra::project(saduria, WGS84)
density_saduria <- terra::extract(saduria_latlon, pred_grid %>% dplyr::select(lon, lat), method = "bilinear")
pred_grid$density_saduria <- density_saduria$FWBiomassm_raster_19812019presHighweightcor_no0_newZi
plot_map_fc +
geom_raster(data = pred_grid, aes(X*1000, Y*1000, fill = density_saduria)) +
geom_point(data = dd_plot, aes(X*1000, Y*1000, size = sample_size), alpha = 0.5, shape = 21, color = "pink") +
scale_size(range = c(.01, 2), name = "# stomachs") +
scale_fill_viridis() +
theme_sleek(base_size = 6) +
facet_wrap(~year, ncol = 5)Warning: Removed 16443 rows containing missing values (`geom_raster()`).
Warning: Removed 2 rows containing missing values (`geom_point()`).
ggsave(paste0(home, "/figures/supp/index_ci_saduria_map.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 16443 rows containing missing values (`geom_raster()`).
Removed 2 rows containing missing values (`geom_point()`).
# What is the average saduria density at the location of the hauls over year?
df_sad <- df
density_saduria <- terra::extract(saduria_latlon, df_sad %>% dplyr::select(lon, lat), method = "bilinear")
df_sad$density_saduria <- density_saduria$FWBiomassm_raster_19812019presHighweightcor_no0_newZi
df_sad <- df_sad |>
group_by(year) |>
summarise(mean_saduria = mean(density_saduria))group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped
p1 <- ggplot(df_sad, aes(year, mean_saduria)) +
geom_point() +
geom_line() +
NULL
p2 <- ggplot(year_cond |> filter(prey == "Saduria"), aes(year, exp(est))) +
geom_line() +
geom_point()filter: removed 58 rows (67%), 29 rows remaining
p1 / p2year_cond2 <- year_cond |> filter(year %in% c(df_sad$year) & prey == "Saduria")filter: removed 65 rows (75%), 22 rows remaining
plot(df_sad$mean_saduria ~ exp(year_cond2$est))# Somewhat good correlation but not perfect. This is not a direct test though, because it assumes they eat saduria in proportion...
# Lastly, calculate the index for the rectangles that we have samples in each year...
# Refit without extra time, because when we filter year_rectangles that are only in data, we don't have all years and cant fit to all years...
fit_sad_m2_extra4 <- sdmTMB(
FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
data = df,
mesh = mesh,
time = "year",
time_varying = ~1,
spatial = "on",
spatiotemporal = "off",
family = tweedie(link = "log"))Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`
summary(fit_sad_m2_extra4)Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
coef.est coef.se
depth_sc -0.01 0.15
spred_length_sc 0.19 0.05
Smooth terms:
Std. Dev.
sds(pred_length_sc) 0
Random intercepts:
Std. Dev.
month_f 0.95
Time-varying parameters:
coef.est coef.se
(Intercept)-1993 -9.19 1.08
(Intercept)-1994 -9.22 1.08
(Intercept)-1995 -8.63 1.01
(Intercept)-1996 -8.61 1.00
(Intercept)-1997 -8.61 1.01
(Intercept)-1998 -9.71 1.11
(Intercept)-1999 -8.58 0.99
(Intercept)-2000 -8.16 0.98
(Intercept)-2001 -7.53 0.98
(Intercept)-2002 -9.07 1.00
(Intercept)-2003 -8.43 0.98
(Intercept)-2004 -8.44 0.98
(Intercept)-2005 -8.24 0.98
(Intercept)-2006 -8.01 0.97
(Intercept)-2007 -8.80 0.98
(Intercept)-2008 -8.67 0.98
(Intercept)-2009 -8.61 0.99
(Intercept)-2012 -7.96 0.99
(Intercept)-2013 -8.48 0.99
(Intercept)-2014 -8.22 0.98
(Intercept)-2018 -7.54 0.95
(Intercept)-2021 -7.14 0.94
Dispersion parameter: 0.73
Tweedie p: 1.54
Matérn range: 80.79
Spatial SD: 2.99
ML criterion at convergence: -782.798
See ?tidy.sdmTMB to extract these values as a data frame.
**Possible issues detected! Check output of sanity().**
df_saduria <- df |> mutate(year_rect_id = paste(year, ices_rect, sep = "."))mutate: new variable 'year_rect_id' (character) with 159 unique values and 0% NA
pred_grid4 <- pred_grid |>
mutate(year_rect_id = paste(year, ices_rect, sep = ".")) |>
filter(year_rect_id %in% c(df_saduria$year_rect_id))mutate: new variable 'year_rect_id' (character) with 1,769 unique values and 0% NA
filter: removed 383,605 rows (90%), 41,941 rows remaining
plot_map_fc +
geom_raster(data = pred_grid4, aes(X*1000, Y*1000, fill = ices_rect)) +
facet_wrap(~year) +
geom_point(data = df, aes(X*1000, Y*1000))Warning: Removed 1249 rows containing missing values (`geom_raster()`).
Warning: Removed 39 rows containing missing values (`geom_point()`).
Warning: colourbar guide needs continuous scales.
pred_sad4 <- predict(fit_sad_m2_extra4, newdata = pred_grid4, return_tmb_object = TRUE)
ncells4 <- filter(pred_grid4, year == max(pred_grid4$year)) |> nrow()filter: removed 33,599 rows (80%), 8,342 rows remaining
index_sad4 <- get_index(pred_sad4, area = rep(1/ncells4, nrow(pred_sad4$data)), bias_correct = TRUE)
manual_index <- pred_sad4$data |>
group_by(year) |>
summarise(index_manual = mean(exp(est)))group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped
index_sad4 |>
ggplot(aes(year, est)) +
geom_point() +
geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = manual_index, aes(year, index_manual, color = "Manual index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = year_cond |> filter(prey == "Saduria"), aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_ribbon(data = year_cond |> filter(prey == "Saduria"), aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
coord_cartesian(ylim = c(0, 0.015)) +
theme(legend.position = "bottom")filter (grouped): removed 44 rows (67%), 22 rows remaining
filter: removed 58 rows (67%), 29 rows remaining
filter: removed 58 rows (67%), 29 rows remaining
# The only thing that could change the conditional from the index is depth I guess? And size
t <- pred_grid4 |>
group_by(year) |>
summarise(mean_depth_sc = mean(depth_sc))group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped
df |>
group_by(year) |>
summarise(mean = mean(pred_length_sc)) |>
ggplot(aes(year, mean)) +
geom_point()group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped
# Make a new conditional with values closer to the prediction grid
nd_year2 <- data.frame(year = unique(pred_grid4$year),
depth_sc = t$mean_depth_sc,
pred_length_sc = 0,
month_f = as.factor(11))
year_pred_sad2 <- predict(fit_sad_m2_extra4, newdata = nd_year2, re_form = NA, re_form_iid = NULL, se_fit = TRUE)
index_sad4 |>
ggplot(aes(year, est)) +
geom_point() +
geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = manual_index, aes(year, index_manual, color = "Manual index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_point(data = year_pred_sad2, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_ribbon(data = year_pred_sad2, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme(legend.position = "bottom")filter (grouped): removed 44 rows (67%), 22 rows remaining
index_sad4 |>
ggplot(aes(year, est)) +
geom_line() +
geom_line(data = year_pred_sad2, aes(year, exp(est), color = "Conditional Year")) +
geom_line(data = year_pred_sad, aes(year, exp(est), color = "Conditional Year fixed")) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme(legend.position = "bottom")# Ok, now I can get the index and the conditional effect match, i.e., by setting depths to match. Can I get conditional to match data by setting depth, month and size to match?
tt <- df |>
group_by(month_f, year) |>
summarise(n = n()) |>
group_by(year) |>
filter(n == max(n)) |>
arrange(year)group_by: 2 grouping variables (month_f, year)
summarise: now 43 rows and 3 columns, one group variable remaining (month_f)
group_by: one grouping variable (year)
filter (grouped): removed 21 rows (49%), 22 rows remaining
t <- df |>
group_by(year) |>
summarise(mean_depth_sc = mean(depth_sc),
mean_pred_length_sc = mean(pred_length_sc))group_by: one grouping variable (year)
summarise: now 22 rows and 3 columns, ungrouped
nd_year3 <- data.frame(year = unique(df$year),
depth_sc = t$mean_depth_sc,
pred_length_sc = t$mean_pred_length_sc,
month_f = tt$month_f) # Most common month!
year_pred_sad3 <- predict(fit_sad_m2_extra4, newdata = nd_year3, re_form = NA, re_form_iid = NULL, se_fit = TRUE)
ggplot() +
geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_line(data = year_pred_sad3, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
geom_ribbon(data = year_pred_sad3, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme(legend.position = "bottom")filter (grouped): removed 44 rows (67%), 22 rows remaining
Warning in geom_line(data = year_pred_sad3, aes(year, exp(est), color =
"Conditional Year"), : Ignoring unknown parameters: `shape`
# Is the mean higher because I have tails that the model doesn't capture?
df_sum2 <- df |>
filter(FR_sad < quantile(FR_sad, prob = 0.95)) |>
group_by(year) |>
summarise(Saduria = mean(FR_sad))filter: removed 622 rows (5%), 11,803 rows remaining
group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped
ggplot() +
geom_line(data = year_pred_sad3, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE) +
geom_ribbon(data = year_pred_sad3, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
fill = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE, color = NA) +
geom_line(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE) +
geom_line(data = df_sum2, aes(year, Saduria, color = "Data trimmed"), alpha = 0.6, inherit.aes = FALSE) +
labs(x = "Year", y = "Feeding ratio", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme(legend.position = "bottom")filter (grouped): removed 44 rows (67%), 22 rows remaining
# YEEEESS FINALLY!!!
# TODO: still I want to estimate the effect of index prediction vs conditional... Plot spatial predictions
spatial_preds <- bind_rows(pred_spr$data |> mutate(prey = "Sprat"),
pred_her$data |> mutate(prey = "Herring"),
pred_sad$data |> mutate(prey = "Saduria"))mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
# Plot spatial random effect
plot_map +
geom_raster(data = spatial_preds |> filter(year == 2000),
aes(X*1000, Y*1000, fill = omega_s)) +
scale_fill_gradient2(name = "Spatial random field") +
facet_wrap(~prey) +
geom_sf(size = 0.1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.1, 0.7)) +
theme(legend.key.height = unit(0.6, "line"),
legend.key.width = unit(0.2, "line"))filter: removed 1,232,616 rows (97%), 44,022 rows remaining
Warning: Removed 1701 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/diet_omega_s.pdf"), width = 17, height = 6, units = "cm")Warning: Removed 1701 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Sprat"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR sprat") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 851,092 rows (67%), 425,546 rows remaining
Warning: Removed 16443 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/sprat_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 16443 rows containing missing values (`geom_raster()`).
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Herring"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR herring") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 851,092 rows (67%), 425,546 rows remaining
Warning: Removed 16443 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/herring_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 16443 rows containing missing values (`geom_raster()`).
plot_map_fc +
geom_raster(data = spatial_preds |> filter(prey == "Saduria"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "FR saduria") +
facet_wrap(~year) +
geom_sf(size = 0.1)filter: removed 851,092 rows (67%), 425,546 rows remaining
Warning: Removed 16443 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/supp/saduria_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")Warning: Removed 16443 rows containing missing values (`geom_raster()`).
# Plot spatiotemporal predictions for a year and all species
plot_map +
geom_raster(data = spatial_preds|> filter(year == "2000"), aes(X*1000, Y*1000, fill = exp(est))) +
scale_fill_viridis(trans = "sqrt", name = "Per cap. predation") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria"))) +
geom_sf(size = 0.1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.08, 0.7),
legend.key.height = unit(0.6, "line"),
legend.key.width = unit(0.2, "line"))filter: removed 1,232,616 rows (97%), 44,022 rows remaining
Warning: Removed 1701 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/spatial_prey_prediction_2000.pdf"), width = 17, height = 6, units = "cm")Warning: Removed 1701 rows containing missing values (`geom_raster()`).
Plot spatial overlap vs and populationa and per capita feeding
# Join feeding ratio indices (per capita) and overlap
cod_spr <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
dplyr::select(year, cod_spr_ovr_tot) |>
rename(overlap = cod_spr_ovr_tot) |>
mutate(prey = "Sprat")Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_her <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
dplyr::select(year, cod_her_ovr_tot) |>
rename(overlap = cod_her_ovr_tot) |>
mutate(prey = "Herring")Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
cod_ben <- read_csv(paste0(home, "/output/cod_ben_sum_ovrlap.csv")) |>
dplyr::select(year, cod_sad_ovr_tot) |>
rename(overlap = cod_sad_ovr_tot) |>
mutate(prey = "Saduria")Rows: 54 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): year, quarter, cod_sad_ovr_tot
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)
mutate: new variable 'prey' (character) with one unique value and 0% NA
overlap <- bind_rows(cod_spr, cod_her, cod_ben)
index_ovr <- index |>
left_join(overlap, by = c("year", "prey")) |>
drop_na()left_join: added one column (overlap)
> rows only in x 6
> rows only in y ( 0)
> matched rows 108 (includes duplicates)
> =====
> rows total 114
drop_na: removed 6 rows (5%), 108 rows remaining
# Join in cod biomass data to calculate snapshot predation
cod_pred <- read_csv(paste0(home, "/output/pred_cod.csv")) |>
filter(quarter == 4)Rows: 943428 Columns: 38
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (3): substrate, month_year, ices_rect
dbl (35): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 471,714 rows (50%), 471,714 rows remaining
# Left join cod biomass density onto the grid of diet predictions
spatial_preds_dens <- left_join(spatial_preds,
cod_pred |> dplyr::select(est, X, Y, year) |> rename(est_codbiom = est))rename: renamed one variable (est_codbiom)
Joining with `by = join_by(X, Y, year)`left_join: added one column (est_codbiom)
> rows only in x 6,177
> rows only in y ( 48,227)
> matched rows 1,270,461
> ===========
> rows total 1,276,638
# Multiply local feeding ratio with biomass density of cod
spatial_preds_dens <- spatial_preds_dens |>
drop_na(est_codbiom) |>
mutate(pred = exp(est) * est_codbiom)drop_na: removed 6,177 rows (<1%), 1,270,461 rows remaining
mutate: new variable 'pred' (double) with 1,270,461 unique values and 0% NA
# Plot map!
plot_map +
geom_raster(data = spatial_preds_dens |> filter(year == "2000"), aes(X*1000, Y*1000, fill = pred)) +
scale_fill_viridis(trans = "sqrt", name = "Pop. predation") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria"))) +
geom_sf(size = 0.1) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.08, 0.7),
legend.key.height = unit(0.6, "line"),
legend.key.width = unit(0.2, "line"))filter: removed 1,226,652 rows (97%), 43,809 rows remaining
Warning: Removed 1698 rows containing missing values (`geom_raster()`).
ggsave(paste0(home, "/figures/spatial_pop_predation_prediction_2000.pdf"), width = 17, height = 6, units = "cm")Warning: Removed 1698 rows containing missing values (`geom_raster()`).
# index? Just take the average
spatial_preds_dens |>
group_by(year, prey) |>
summarise(mean_pop_pred = mean(pred)) |>
ggplot(aes(year, mean_pop_pred)) +
geom_point(alpha = 0.7) +
stat_smooth(method = "gam", formula = y~s(x, k=3), color = "steelblue") +
facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
labs(x = "Year", y = "Population predation", color = "") +
scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
theme_sleek(base_size = 9) +
theme(legend.position = c(0.9, 0.9))group_by: 2 grouping variables (year, prey)
summarise: now 87 rows and 3 columns, one group variable remaining (year)
ggsave(paste0(home, "/figures/pop_pred_index.pdf"), width = 17, height = 7, units = "cm")
# Summarize across years
spatial_preds_dens_sum <- spatial_preds_dens |>
group_by(prey, year) |>
summarise(tot_pred = sum(pred))group_by: 2 grouping variables (prey, year)
summarise: now 87 rows and 3 columns, one group variable remaining (prey)
# Left_join with index data
index_ovr <- index_ovr |>
left_join(spatial_preds_dens_sum, by = c("year", "prey")) |>
drop_na()left_join: added one column (tot_pred)
> rows only in x 0
> rows only in y ( 6)
> matched rows 108
> =====
> rows total 108
drop_na: no rows removed
# Check some regressions...
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Herring")))
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Saduria")))
#
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Herring")))
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Saduria")))
# Calculate correlations between FR, predation and overlap
cor <- plyr::ddply(index_ovr, c("prey"),
summarise,
cor_fr_ovr = round(cor(overlap, est), 2),
cor_pred_ovr = round(cor(overlap, tot_pred), 2)) |>
pivot_longer(c("cor_fr_ovr", "cor_pred_ovr")) |>
mutate(name = ifelse(name == "cor_fr_ovr", "Per capita", "Population"))summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
pivot_longer: reorganized (cor_fr_ovr, cor_pred_ovr) into (name, value) [was 3x3, now 6x3]
mutate: changed 6 values (100%) of 'name' (0 new NA)
p <- index_ovr |>
rename(fr = est) |>
pivot_longer(c("fr", "tot_pred")) |>
mutate(id = paste(name, prey, sep = ";")) %>%
split(.$id) |>
purrr::map(~lm(value ~ overlap, data = .x)) |>
purrr::map_df(broom::tidy, .id = 'id') |>
filter(term == 'overlap') |>
separate(id, sep = ";", into = c("name", "prey"), remove = FALSE) |>
mutate(name = ifelse(name == "fr", "Per capita", "Population"))rename: renamed one variable (fr)
pivot_longer: reorganized (fr, tot_pred) into (name, value) [was 108x10, now 216x10]
mutate: new variable 'id' (character) with 6 unique values and 0% NA
filter: removed 6 rows (50%), 6 rows remaining
mutate: changed 6 values (100%) of 'name' (0 new NA)
p# A tibble: 6 × 8
id name prey term estimate std.error statistic p.value
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fr;Herring Per capita Herring over… 3.06e-4 3.11e-3 0.0982 0.923
2 fr;Saduria Per capita Saduria over… -2.48e-4 4.88e-4 -0.509 0.613
3 fr;Sprat Per capita Sprat over… 9.58e-3 4.88e-3 1.96 0.0610
4 tot_pred;Herring Population Herring over… -2.83e+1 1.39e+4 -0.00203 0.998
5 tot_pred;Saduria Population Saduria over… 1.31e+3 5.64e+2 2.32 0.0243
6 tot_pred;Sprat Population Sprat over… 2.62e+4 1.94e+4 1.35 0.188
# Plot!
index_ovr |>
rename('Per capita' = est,
'Population' = tot_pred) |>
pivot_longer(c('Per capita', 'Population')) |>
ggplot(aes(overlap, value)) +
geom_point(alpha = 0.5) +
ggh4x::facet_grid2(name ~ prey, scales = "free", independent = "y") +
# geom_text(data = cor,
# aes(label = paste("r=", value, sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
# inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
geom_text(data = p,
aes(label = paste("p=", round(p.value, digits = 3), sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
labs(y = "Predation", x = "Spatial overlap", color = "Year") +
theme_sleek(base_size = 9) +
theme(legend.position = "bottom", aspect.ratio = 1) +
NULLrename: renamed 2 variables (Per capita, Population)
pivot_longer: reorganized (Per capita, Population) into (name, value) [was 108x10, now 216x10]
# Summarize the differences here, make fake facet grid free axis
ggsave(paste0(home, "/figures/fr_overlap_cor.pdf"), width = 17, height = 11, units = "cm")